The cosmological behavior of Bekenstein's modified theory of gravity 
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I. INTRODUCTION 

Bekenstein has proposed a covariant and relativistic 
theory of gravity [l| which is to meant to provide for an 
alternative to dark matter. The theory is called TeVeS 
(which stands for Tensor- Vector-Scalar) and can have 
the Bekenstein-Milgrom Aquadratic Lagrangian 2] non- 
relativistic limit, whichprovides for Milgrom's Modified 
Newtonian Dynamics [j] (MOND) as an explanation of 
galactic rotation curves. 

MOND has had a number of successes at fitting data 
on galactic scales, less so on the scale of clusters. As 
shown by Sanders Q, a universe governed by MOND 
still requires a small fraction of massive neutrinos (corre- 
sponding to a neutrinos mass of approximately 2 eV) to 
account for the dynamics of clusters. More recently, the 
data obtained in Clowe et al 0] indicates that some extra 
mass would be needed in the "Bullet" cluster to account 
for the mismatch between X-ray and weak lensing recon- 
structions of its mass profile. Again, this can be shown, 
in the context of MOND, to be resolved with a small 
fraction of massive neutrinos, as claimed by Sanders i. 
However none of these results need to be valid in TeVeS: 
the vector field which was shown to be extremely impor- 
tant for structure formation Q, is also expected to be 
important on the scale of clusters. 

MOND is the nonrelativistic limit, in the regime of 
small acceleration, of TeVeS. This is exactly akin to New- 
tonian gravity being the nonrelativistic limit of General 
Relativity. And in the same way, it is incorrect to in- 
fer all properties about the relativistic theory from the 
nonrelativistic regime. The most notable case is in gravi- 
tational lensing: the Newtonian limit supplies the wrong 
answer if one wishes to study general relativity. Hence 
results obtained in the MONDian regime do not neces- 
sarily extend to TeVeS and in many situations one must 
explore the fully relativistic theory. One regime where 
this is imperative is the large scale dynamics of the Uni- 
verse. In this paper we try to infer general properties of 
a universe governed by TeVeS. 

TeVeS incorporates a dynamical scalar field 4>, two 



metric tensors and a unit-timelike vector field A". In 
the case of a homogeneous and isotropic space-time, 
the TeVeS field equations reduce to a form similar to 
Friedman-Robertson- Walker (FRW) form. There have 
been several studies [H, 0, H, S [l^ of the resulting back- 
ground cosmology with different choices of the TeVeS free 
function. In this paper we extend the preliminary anal- 
ysis of TeVeS cosmology to more general settings. 

The paper is organized as follows. In Sec. 2 we give the 
general TeVeS equations, their applications to cosmology 
and the notation we will use throughout the paper. In 
Sec. 3 we will present our numerical results and we will 
discuss their behavior with some analytical insights. We 
conclude in Sec. 4. 



II. FUNDAMENTALS OF TEVES COSMOLOGY 
A. TeVeS action 

Let us first describe the basics of TeVeS theory. We 
use the conventions of [llj . 

In TeVeS we have two tensor fields which are the two 
metrics g^u and g^^i/ mentioned above, a scalar field and 
a vector field A°'. They can be linked through 

9^.^e-^''{g^,+A^A,)-e^^A^A, (1) 



The total action of the system is Sg- 
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Sm where Sg is the traditional Einstein-Hilbert action 
built with the Einstein frame metric , Sg is the action 
for (f) and /i, Sy is the action for the vector field and Sm 
denotes the action of ordinary matter. Explicitly we have 
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where G is the bare gravitational constant, Kb is a 
dimensionless constant, A is a Lagrange multiplier im- 
posing a unit-timelike constraint on the vector field, 
Fap — 2V[q^^] and /i is a nondynamical field with a free 
function V{^). The matter Lagrangian L,„ depends only 
on the matter frame metric g^j,i, and a generic collection 
of matter fields ■ 

Although the theory was designed to provide for a rel- 
ativistic theory of MOND, taken at face value this need 
not be the case. It all lies within the choice of the func- 
tion V{ijl). The requirement that both exact MOND and 
exact Newtonian limits exist within the theory puts some 
constraint on the possible forms of V{fi). In particular 
these two requirements fix the derivative of V{ijl) to be 
of the form 
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where £b is a scale, a- dimensionless constant and /(/i) 
is still an arbitrary function, not related to the MOND or 
Newtonian limits, whose only constraint is that /(O) = /q 
is nonvanishing. The in the numerator is what is pre- 
cisely required to have an exact MOND limit (which is 
reached as /i — > 0), whereas the factor 1 ~ fi/fiQ is what 
is required to have an exact Newtonian limit (which is 
reached as /x fj,o). Bekenstein chooses the simplest 
case, namely m — 1, which we shall also do in this work. 
The measured Newton's constant Gn and MOND accel- 
eration oo are given irrespectively of the form of /(/Lt) or 
the power m, as 
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generalize this choice as 



(5) 



where /t = /j,//j,o, /ia is a constant and c„ a set of addi- 
tional constants. We also let n run over negative values. 
The Bekenstein model is recovered as /i^ = 2, C2 = 1 and 

Cn^2 = 0. 

Different functions have also been proposed especially 
by and and slight modifications to Bekenstein's 
toy model have been given by [To| . 



B. Homogeneous and isotropic cosmology 

For the purpose of this paper we will restrict ourselves 
to homogeneous and isotropic spacetimes (see [ll| for 
a detailed investigation of the inhomogeneous case). We 
adopt the Friedmann-Lemaitre- Robertson- Walker metric 
with the scale factor a for the matter frame metric and 
b = ae'* for the Einstein frame metric. We choose a 
coordinate system such that the components of the vec- 
tor field take the form A" = (1,0,0,0). The modified 
Friedmann-equation becomes 



3H^^87rGe-^^ {p^ + Px) 



(6) 



where H ~ h/h is the Einstein-frame Hubble parameter 
with an overdot means a derivative with respect to the 
coordinate time chosen and where we have defined a field 
density as 
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and 



The ordinary fluid energy density px evolves as usual as 
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where 0c is the cosmological value of (j) when the qua- 
sistatic system in question broke away from the cosmo- 
logical expansion. 

What remains is to choose the (still free) function f{p). 
It turns out that for cosmological models we must have 
V' > Q. Therefore in the m — \ case that we are con- 
sidering, this cannot happen in the same branch as for 
a quasistatic system, i.e. we need p, > p,o for cosmol- 
ogy. (There is a very clever alternative considered in 
[i^l with its cosmology studied in We do not follow 
that approach here). In the case that there is an single 
extremum in V' for p > po one has a further choice im- 
posed by the requirement that V' be single valued in the 
branch considered : either use the branch from p — po 
up to the position of the extremum or the branch from 
the extremum to infinity. Bekenstein makes the second 
choice, which is what we will also do in this work. 

With the above considerations in mind, Bekenstein 
makes the choice f{p) = {p/po — 2)^. We will therefore 



Px = -3-(l + w)px 
a 



(8) 



where w is the equation of state parameter of the fluid. 
Relative densities fli can be defined as usual as 



Pi + PcP 



(9) 



Variation of the action with respect to the field p gives 
back the constraint 



(10) 



which fixes p as a function of (j). 

Finally the scalar field evolves according to a system 
of two first order equations given by 
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and 



t + 3HT^8TTGe-^'''px{l + iw). (12) 
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III. NUMERICAL RESULTS 



A. Bekenstein's toy model 

Let us first revise the results presented in i.e. that 
(j) acts as a tracker field in the radiation, matter and A 
eras. 

In this special case of the Bekcnstein toy model, the 
function V can be explicitly inverted analytically by 
setting q = 647r?|0V(3Aio), r = g^g7l2 + 64/81 s = 
(27q2 + 128)/54,pi = (r+s)i/3,p2 = V9pi + 12 + 16/pi, 
P3 — (9pi — 24pi + 16) p2 — 5iqpi and one obtains an ex- 
plicit expression for the nondynaniical field /i: 



1 

^=6 
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(13) 



In figdiwe plot il^ as a function of log(a) for the dif- 
ferent epochs in the Universe history: panel a) radiation 
era, panel 6) matter era and panel c) the A era. As can be 
seen in the figure, for several different initial conditions, 
(j) has an attractor in the different epochs, which leads to 
the tracking behavior. 
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FIG. 1: f2<^ as a function of log(a). The units for the scale 
factor are arbitrary. has an attractor which correspond 
to tracking the three constituents independently of the initial 
conditions. We illustrate this here, taking /io = 4 and plotting 
three different histories in: a) radiation era, fe) matter era and 
c) the A era. 

In panel (a) of figE] we plot the fractional energy den- 
sities fii as a function of log(a). As can be seen from the 
figure, (/) synchronizes its energy density with the dom- 
inant component of the Universe. As mentioned in Q, 
during the radiation era 17^ reaches 2^ a-nd during the 
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FIG. 2: Panel a) represents Bekenstein's toy model, panel and 
6) is a model with C2, c-i and C-2 as non vanishing. The long 
dashed line represents the dot-dashed one is r2A(i), the 

solid line is and the short dashed line is for JIm (i). 



matter and A eras the limit is — ^From panel (a) 
of figE] it is clear that the overall behavior of (j) is simply 
a sequence of the three tracking epochs. This behavior 
is similar to other general cosmological theories involving 
tracker fields llSl. 



B. Generalized function 

The behavior of (f) depends on the form of the function 
/(/i) chosen. In this section we investigate the main fea- 
tures which result from a function given by equation ^ . 

In FigO we plot the evolution of fi^ in the presence 
of nonrelativistic matter. With only C2 as nonzero co- 
efficient, we have an apparently stable attractor. If we 
choose in ([5]) a nonzero Ci (while keeping the nonzero C2) 
we find that the tracking is still possible but less stable 
than the one with C2 alone. Hence, we can have longer 
lived tracking by decreasing ci . Adding cq with suitably 
chosen coefficients (see subsection on analytical results 
below) does not change this behavior as can be seen in 
FigEl The curves corresponding to the two cases just 
discussed essentially overlap, and depart from tracking 
at log a = —4. It seems that, in general, mixing cases 
of different n together will lead to the loss of tracking 
behavior, with at best having temporary trackers. While 
this turns out to be the most frequent behavior there are 
counter examples : a mixed case of n = 1 and n = 3 
with fia = 5/7 and C1/C3 = —4/49 also exhibits perfectly 
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FIG. 3: vs log a for the generalized function given by equa- 
tion ([5]). The solid line corresponds to C2 7^ 0. The dotted 
and dotted-short-dashed curves mixed together correspond to 
(c2,ci) and (c2,Ci,co) as the only non vanishing coefficients. 
The long-dashed curve corresponds to (c2 7^ 0,c_i 7^ 0) and 
the dotted-long-dashed one to (c2 7^ 0,c_i 7^ 0, c_2 7^ 0). 
Finally the short-dashed curve on the left corresponds to 
C-2 + 0. 

stable tracker behavior (see fig. 2]). As we will see in 
the next section on analytical treatment, the existence of 
tracking behavior or not, can be understood in terms of 
a generic set of rules. 

Tracking behavior disappears altogether when we add 
a negative power in (/i — /Iq). We can take either c_i 7^ 
or C-\ 7^ and c_2 7^ and the effect will essentially 
be the same. The curves corresponding to these cases 
almost overlap in Fig[3](by chance), and fi^ rapidly tends 
to 1 at log a = —8. In both of these cases, the system 
exits tracking quite quickly. The initial stage tracking 
is a result of the presence of C2; if we set C2 = and 
keep c_2 7^ 0, for example, the leftmost curve of Fig[3] 
indicates that tracking is impossible. Higher powers of 
(Ji — Ha) can suppress this instability temporarily. 

What if we add more negative powers of n? The lower 
panel of Fig[5]shows the evolution of the relative densities 
fii in the case where only c_3 7^ and where only a 
cosmological term is present. 

The middle panel of the same figure show that T/H is 
a constant during tracking while (j) tends to a constant in 
the infinite future. The upper panel shows that H also 
tends to a constant. In the case where only radiation 
and/or pressureless matter is present, we find that evo- 
lution has an early unstable tracking behavior just like 
the A case, but eventually evolution stops in finite time. 
Adding a cosmological term can remedy the situation by 



FIG. 4: A mixed case with n = 1 and n = 3 which exhibits a 
perfectly stable tracker, in the radiation era. 
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FIG. 5: A case with only c-3 7^ and a cosmological constant. 
Upper panel : Evolution of ln(///J/oo) withlnb (see section on 
the analytical treatment of the n = — 3 case for the definition 
of Hoc) ■ Middle panel : Evolution of <l)/4'oo (sohd) and T/H 
(dotted) with In b (see section on the analytical treatment of 
the n — —3 case for the definition of (l)oa)- Lower panel : 
Evolution of the relative densities fi^ (dotted) and JIa (solid) 
showing an initial temporary tracking phase followed by A- 
domination. 
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creating a second tracking era. 



/Za(l + e) with e <C 1. The impact of different terms in 
the function series is investigated bellow. 
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FIG. 6: A case with only c_4 7^ 0. Upper panel : Evolution 
of In /i (solid) with In b. Middle panel : Evolution of (p (solid) 
and T/H (dotted) with In 6. Lower panel : Evolution of rel- 
ative densities of radiation (solid) and scalar field (dotted) 
versus In 6 with only a radiation fluid present, showing a late 
radiation domination. 

Considering even more negative terms, e.g. only c_4 ^ 
as in Fig [6] gives yet different behaviors. An early track- 
ing era eventually conies to an end followed by a fluid 
dominated era (fix = 1). During both eras we find that 
the field F always evolves linearly with log t (if the back- 
groimd fluid is radiation). 

Our numerical study of the dynamics of this system 
picks out four types of behavior: perfect tracking, track- 
ing followed by domination of fi^, no tracking and even- 
tual singular behavior and finally initial tracking and 
eventual fluid domination. Given that we have consid- 
ered a reasonable generalization of the function proposed 
by Bekenstein, we can see that these types of behavior 
will be widespread in TeVeS. We can now proceed to un- 
derstand what lies behind these different regimes. 



IV. ANALYTICAL RESULTS 

In this section we investigate the existence and stability 
of tracker solutions and develop approximate analytical 
solutions to the evolution of the scalar field for the various 
choices of function and tracking regimes. 

When does the field track? The numerical studies of 
the previous section showed that tracking occurs as V 
tends to its zero point. It therefore makes sense to study 
solutions when fi is close to this point, i.e. we let (i = 



A. Case n > 1 

In this case the function V and its derivative V' are 
given by 



Mo 
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and 



m— 1 

+ (l-Aia)"ln|A-l| 



(n+ l)(n + 2) 

-m 

—(A -Mar 



(M - Ma)' 



(14) 



(15) 



respectively, where we put a subscript (n) to denote the 
case considered, and where /„ is an integration constant 
(for the case n). 

Clearly Vj'„)(Ma) = for any n. In the case of V(„) 
however, it depends on the choice of the integration con- 
stant /„. We will discuss the relevance of this constant 
further below. 

Let us choose for the moment the constant /„ such that 
V(„)(/ia) ~ 0. We then expand the relevant quantities as 
a Taylor series about the point /I = /^q. We have 



and 



y ^ MoMaK) 
71+1 



V = Voe" 



(16) 



(17) 



where % = (11(11+"- /{l6T:e%{^ia - 1)). The field density 
is therefore given by 
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to lowest order in e. 

Using the constraint equation which links fi with (/), 
namely PU)) . and taking the square root, we obtain 



(19) 



where /3 = ±1 is an artifact of taking the square root. 
The Taylor expansion of the auxiliary scalar field F is 




(20) 
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Since we are interested on finding possible tracking be- 
havior, motivated by the numerical results, it is reason- 
able to make the ansatz — Apx with A a constant. 



givmg 



Using the Friedmann equation we get 



Px 



(21) 



(22) 



Eliminating e" from the auxiliary scalar field F we get 



F = -2f3y^8TTGAna 



(23) 



Using (|12p along with the derivative of (|23p we find 



A 



(l + 3w)2 



1 + A 3paPoil~wy 



(24) 



More specifically in the three special eras considered 
above, namely radiation, matter and A eras, we get 



in matter and A eras, while fJw 



2 for the Bekenstein 



in the radiation era. Taking 

toy model we obtain the results previously found in ^ 

Now we specify the constant /3 which tells whether (p 
is negative or positive. Changing the time variable to 
In a and using (|22p and ITOl) , one finds that = 4>i is 
approximately a constant during tracking, given by 



1-K3w 



MaAio(l - w)l3 - (1 + iw) 



(25) 



This results to 6 — 



bi Ina. Introducing the above 



2i/ii+3(l + m) 

in (I2ip we find e cx a " where po is the fluid 
density at a = 1. Consequently e will be a decreasing 
function going to asymptotically, and e ^ if and only 
if tracking is stable, that's to say if 2(/)i -f 3 (1 + w) > 0. 
Using (|25p this condition is equivalent to 



-{l + 2,wf + ?,{l~w^) PaPol3 
MaMo (1 - w) /3 - (H- iw) 



> 



(26) 



The above inequality is satisfied if and only if either P — 1 
and w > —1/3 or /? = — 1 and w < —1/3. Therefore /3 
is completely fixed by the equation of state of the back- 
ground fluid. In particular we see that (j> has to change 
sign when passing from the matter to A era which results 
to p0 going momentarily to zero [sj. 

Now let us discuss the role of the integration constant 
/„. Its sole impact would be to add a constant in the 
Friedmann equation. What can motivate such a con- 
stant? For exact MOND limit we need y as /i ^ 0. 
It is in general impossible to have both V{Q) — and 
Vi^J■a) = 0, except in a very special and unique case : the 
Bekenstein toy model, i.e. n = 2 and pa = 2. In every 
other case (not a mixed case with many c„ nonzero, which 
will be consider in another section), setting V{fia) = 



would spoil the exact MOND limit. The impact of such 
a constant in the evolution equations however, is what 
precisely destroys the perfect tracker. In other words as 
long as /„ ^ 167rG'e~^'^px, we would have the tracker so- 
lution found above. When eventually px decreases below 
this threshold, which is bound to happen always, track- 
ing stops. In such case there are two possible behaviors 
after tracking stops. If V{pa) < the requirement that 
the LHS of the Friedmann equation is positive, creates a 
pathological situation where fi^ decreases without bound 
below zero while the fluid relative density increases un- 
bounded above one to counterbalance. In the opposite 
case where V{pa) > the universe enters a de-Sitter 
phase in the Einstein frame where b/b tends to a con- 
stant, and fl^ tends rapidly to 1. This second case is 
what stops the temporary trackers found in the previous 
section numerically, in the mixed cases. 

We end this case by noting that the Hubble parameter 
in our approximation will be slightly different from the 
one we obtain in General Relativity but rather if the 
background radiation temperature is T, it will be given 
hy H (X T'y where 7 = 20i |(1 + w). This can be 
potentially important when calculating nucleosynthesis 
constraints, as in the radiation era we would have 7 = 
2(1 + (pi) rather than the canonical value 7 = 2. The 
impact of this kind cosmologies on nucleosynthesis has 
been investigated in [l6l |. 



Case n = 



In this case the function is given by 
1 



16n£ 



^ \^o+ +A + ln|/i-l| ) (27) 



and 



167r£| A - 1 ■ 

The constraint can be inverted analytically to give 



Mo 



4:^/^T > 



(28) 



(29) 



It is therefore clear that \(j>\ is bounded from below : If/)] > 
Po/{V8tt£b)- At this value p. takes its minimum possible 
value given hy p.c = 2. The upshot is that contrary to 
the case where n > 1, we have X^'(/ic) = Mo/(4'''^s) i-^- 
it is not zero. We will show in this subsection that this 
creates a pathological situation. 

First note that since cannot be zero, the two branches 
of positive and negative cf) are disconnected and we can 
change variables from (p to fi. The equivalent of (fT^ is 
then given by 



-p + 6Hp-\ 2 : y = (30) 



A-1 



MS 
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where y = SirGe '^'^px and (i is once again the sign of 0. 

Suppose that the fi — 2 point is reached at some mo- 
ment in time t^. Then V{fi = 2) = ^g(/o + 12)/(167r£|) 
giving 



8TTi% 



2/3(1 + 3w)y/^£i 



M5 



'V2 



(31) 

where the subscript " c" on any variable denotes the value 
of that variable at time tc, and where we have assumed 
that b/b > 0. We consider again two subcases depending 
on the sign of /3(1 + 3u;). 

If /9(1 + 3w) > 0, then jl is always negative (including 
the point jlc) except the very special case where Jq = —12 
and ?/c = 0. Since however we are free to choose the initial 
condition for ji and y independently, we will get that 
generically p,c < 0. Therefore the system of differential 
equations is bound to become ill-defined at some point 

If /3{l+3w) < then there are values of /t (and also /tc) 
which are positive given a large enough y, which means 
that jj, will bounce off the point fx — 2 and start increasing 
again. However as the equations are propagated further 
in time, y decreases even further until once again fl < 
and jj, will decrease back to /i = 2. This time though 
jlc < 0, and the system of differential equations becomes 
once more ill-defined. 



2. Case n ■ 



The equivalent of ^T2\\ is then given by 



2(A-1)(A-Ma) 



2/12 _ 

4/3(1 



3(/1q + l)/i 
I- 3w)^£b VTA 



3Hfi 



1)(A - Ma) 



M0V2 
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(35) 



where y = SnGe '^'^px and /3 is once again the sign of </>. 
Now for /t > we have that 2p? — 2>{pa-\-^)p+^P'a = at 
a value p. ^ p^ ^ \ 3{pa + 1) + y/Qpl - ^'ipa + 9 • In 

other words there the evolution of p reaches a singularity 
at p — pc- To make things worse, since p is negative 
for p > Pc while being positive for p < pc then this 
singularity is always reached! This is true independently 
of the sign of (3 for the same reason as the case n = 0, 
i.e. for positive H, as y decreases below some threshold 
value the term proportional to y can be neglected. 



3. Case n = -2 

Let us now discuss the n = —2 case. 
In this case the function is given by 
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1- Pa P - Pa [Pa - 1)^ 
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In \p - 1 



In \p - pa\ 



(36) 



It turns out that for negative n we need to consider 
each case separately for the first few n's. Lets start with 
the n = — 1 case. 

In this case the function is given by 



^(-1) = 
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In \p — 1 
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and 



^(-1) = 
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The constraint can be inverted analytically to give 



(33) 



{Pa + l)4>^ + \^\J{pa-l?r + 



8^ 
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(34) 



Restricting the allowed range for p to be greater than pa 
(since we have a singularity in V' a,t p — pa) we see that </> 
is once again bounded from below as \4>\ > pq/ {AIbV2tt), 
similarly to the n = case. Since the two branches of 
positive and negative 4> are once again disconnected, we 
can change variables from </) to /t as in the case n = 0. 



and 



V! 



Mo 



(-2) 16^^B (M-1)(M-Ma)2 



(37) 



Contrary to the last case ( n = — 1), is no longer 
bounded from below and can take all possible values. 
However the point </> = can only occur at /i — > cxd, 
which means that once again the two sectors of positive 
and negative (j) are disconnected, since the /i ^ cx3 is a 
singular point (the Friedman equation blows up). We 
can therefore follow the approach we took in the n = — 1 
case and change variables from (j) to p giving 



M = 



2{p-l){p-pa) \^ 
p^ - {3pa + 2)p + ipa I 

4/3(1 -I- 3w)^£b{p - Pa)Vp^ ,^ 
PopV2 



(38) 



where y — SnGe^^'^px and /3 is once again the sign of (p. 
We immediately see that the same problem as the n = — 1 
case arises : there is a singularity in the evolution when 

r 3pa + 



P = Pc 



4- y/ 9pl - ipa + 4 . Since pc > Pa 

always while at the same time p < for p > pc and 
p > for p < Pc , then once again this singularity is 
always reached. 
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4- Case n = —3 

We now turn to the n = —3 case. 
In this case the function is given by 



^(-3) = 



Mo 



/-3 + 



1 



IGniBV ' 2(1 -/i,) (A- /ia)2 

^J,ai^J.a - 2) 1 
(1 - ^a)^ A - /^a 



1 



In 



A-1 



(1 - fla)^ fj. - ^J.a 



and 



^(-3) = 



167r^B (A - 1)(A - Aia)^ 



(39) 



(40) 



Contrary to the n — —1 and n = —2 cases, in this 
case, 2V' + ^V" never goes to zero for fi > ^a, hence 
no singularity of the same type as n = — 1 and n = — 2 
occurs in this case. We consider two hmiting cases of 
behavior. 

The first case is when A approaches /i^- In this hmit 
we get that V' is of order cjp while V is of lower order, 
giving fiV' + V = 2fiQfia<P^- Therefore in this limit the 
equations take the form 



3H^ ^ MoMa^^ + SttGc-^'^px 



and 



(l + 3u;), 

2^io^J■a 



(41) 



(42) 



First we consider the case when the fluid is a cosmo- 
logical constant. Letting e and e to go zero, we get 
that the Hubble parameter tends to a constant H 
Hoo = /ig^/2/(80F^B), while = ^l^^^^ 

and P = 1. Perturbing (f> about (f>ac as = 0oo(l + S) 
we get that = ^^(1 - 2000(5). Differentiating ([44]) 
once under these approximations and eliminating S via 
2mo0oo<5 = Hoc e we get 



{3fia + 1)mo 



e = 



(45) 



This has decaying solutions of the form e = eie "^+* 
£26"^-* where 



1 ± 



3(3/Xa + 1)mo 



(46) 



provided that {S^a + ^)^^o > 8/3. The final stage of this 
particular case is a universe with JIa = 1- 

It turns out that when the fluid is not a cosmological 
constant then, the evolution reaches e — > in finite time. 
This is once again a singularity of a similar form to the 
n = case. The reason we have such a singularity is 
the same as the n = case : one can choose initial con- 
ditions such that (j) — which corresponds to /Lt — > oo, 
at any point in time. However this is inconsistent with 
the equations of motion except in the special case of the 
presence of a cosmological constant. Therefore in this 
case, if the evolution is to be nonsingular in the infinite 
future, one must include a cosmological constant. 



in other words the system behaves as a one with a canon- 
ical scalar field coupled to matter. In this limit we have 

H, which 

H 



tracker solutions such that 
gives VtA, = 



Viw) 



(l-lu)^tOMa 

This tracker however is unsta- 



3(1 — tOj^/^O/Ja 

ble (but can be long-lived depending on the fluid energy 
budget), i.e. as p, becomes larger this limit ceases to be 
valid. 

The second case is when ji oo. We ex- 

pand all functions in powers of e = 1/A, giving 
V' = ^lle^ [1 + (3/ia + l)e] /(167r^|) resulting to low- 
est order in e, = (3^f)e/ {Mb^/^t:), and F = 
-(3^xl^/2 [1 + (3/ia + l)e/2] /(4^b y^). It also turns out 
that although yV' is of 0(e), this order is precisely 
canceled from a term in V , resulting to + V = 
Mo(3/ia + l)e^/(327r£^) to lowest order in e. The evo- 
lution equations become 



Mo 



647r^| 



(43) 



and 

3/Xa + 1 



-3 1 



3/ia + 1 



H- 



4(1 + 3w)(3iBV^V 



M5 



IV2 



^ 
(44) 



5. Case n < —4 

Finally we consider the n < —4 case. 
In this case the function is given by 



Mo 



I-r 



M? 



1 



{n - 1)(1 - /ia) (A - Ma) 



ri-1 



flaif^a - 2) 



1 



n— 3 



1 



(l-Ma)"^^m(A-Ma)' 



In 



M 



(1 - /ia)" A - Ma 



and 



Vf 



(-») 



Mo 



M 



(47) 



(48) 



16^^^ (A -1)(A- Ma)" 

It is easy to show that 2V' + ^V" ^ for A > Ma unless 
A ^ oo. We therefore expect to have smooth evolution 
and furthermore we expect that A will eventually increase 
to infinity. 
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Eventually fi increases to large values where we can 
expand V and V in powers of e = For V' we get 



V 



(49) 



giving 



V2 



(50) 



where /3 is the sign of F. 

It turns out that piV' + y is at least of O(e^) for n = 
—4 or higher for n < —4 while (/> is at least of 0(e^/^) 
for n = —4 or higher for n < —4. We can therefore 
assume that the universe evolves like Einstein gravity 
with H = H and we have the usual Friedmann equation 
: 3i?2 ^ SnGpx 

We consider three separate subcases : w ~ —1, w ^ 1 
and — 1 < w < 1. 

For the w ~ —I case we have that H will tend to a 
constant. The evolution equation for F reads 



F + 3iJF + 6H^ = 



The solution is 



F ^ 2H{e 



-3Ht 



1) 



(51) 



(52) 



and therefore F ^ as required by consistency. Feeding 
back to ([so]) we get that (3 = ~1. 



For the w = I case we have that H = j which means 



that F evolves according to 



The solution in this case is 



6 Fo 

t <3 



(53) 



(54) 



B. Mixed cases 

Mixing cases for n > 1 can have two generic kinds of 
behavior depending on whether V{^a) = or not. 

If V{^a) = then we have perfect trackers. If 
V^Ha) > there will be a temporary tracker until 
V{fJ'a) ~ l6TTGe~^'^px after which the tracker is de- 
stroyed and we get 4> domination. The solution in the 0- 
domination era is de-Sitter space (in the Einstein frame). 

Now the condition of exact MOND means that only a 
very restricted set of function can have V{pa) = 0. For 
the single n cases only the Bekenstein toy model (n — 2 
and /ia = 2) has this feature. By mixing cases however, 
this feature can arise again for very specific choices of fia 
and n's. For example this is impossible for a mixed n = 1 
and n = 2 case, but possible for a mixed n = 1 and n = 3 
case with pa = 2/7 and ci/ca = —4/49. 

While individually n = ~\ and ti — 2 lead to singulari- 
ties, mixing the two can for some choices of c_i and c_2 
lead to 2V' + pV" being nonzero for all fl > fia- This will 
lead to a better behaved evolution although problems can 
still persist as fj, oo just like the n = —3 case. This 
last problem could be removed by mixing in a positive 
power of n. 

Mixing n = with positive n alone cannot remove the 
problem when fi ^ fia^ i-e the evolution will still become 
ill defined at that point. One can still however mix n = 
with both positive and negative powers; the effect of the 
negative power is to drive fl away from the jl = jia point. 

Finally the most general mixed case including both 
positive and negative powers with suitably chosen coeffi- 
cients Ci, will give a function V' which goes to zero not 
ai fl = Pa but at a shifted point fi = p!^. the large p 
behavior will be dominated by the positive powers of n 
while the small p behavior by the negative powers with 
the function still being singular at p = pa- If we restrict 
the evolution for values p > p'^, then we have a situation 
similar to the positive n cases. We therefore expect to 
have the same behavior as if we had a function expanded 
about fl = p'a with positive powers only, i.e. we expect to 
get stable trackers or in the case of a nonzero integration 
constant in V eventual (A-domination. 



where Fq is a constant. Once again F — > as required by 
consistency. Feeding back to ([50]) we get that (3—1. 

The last case is when —\<w<l. In this case F 
evolves as 



t + mv - i(l + Zw)H'^ =0 

which has solution 

^_ 2{l + 'iw) ^^A{l + ?>w) 1 



1 — w 



3(1 -w2) t 



(55) 



(56) 



and therefore F consistently goes to zero. Feeding back 
to (inOl) we get that /3 = 1 for w > -1/3 and (3 = -I for 
w < -1/3. 



V. CONCLUSIONS 

In this paper we have tried to extract general prop- 
erties of the cosmology of TeVeS gravity. We have re- 
stricted ourselves to analyzing the background equations 
for a homogeneous and isotropic universe. We have con- 
sidered a general function for p which we believe encom- 
passes most of the possible regimes one could encounter 
in this theory. 

We find that Bekenstein's choice of V{p) naturally 
leads to stable tracking behavior. This confirms the re- 
sults of Q. Generalizing the Bekenstein function by in- 
cluding arbitrary strictly positive powers of {fi — pa) we 
either get exactly the same type of tracking or tempo- 
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rary tracking followed by (/)-domination. What controls 
this kind of behavior is the value of the integration con- 
stant in V{ii). If V{iiofia) = then the tracking behavior 
is retained to the infinite future, otherwise the tracking 
period is temporary followed by ^-domination. 

We find that the rt = 0, n = — 1 and n = —2 cases 
alone lead to future singularities in the evolution. The 
same happens in the n = —3 case unless the background 
fluid is a cosmological constant. For n < 4 the evolution 
is nonsingular and eventually leads to fluid domination 
with fj, oo. 




2 3 



FIG. 7: Phase plot for the Bekenstein's model. All the curves 
join the tracker solution (57^=0.4) and tracking remains sta- 
ble. 

Finally, the singularities can be avoided by mixing 
cases together. If positive and negative cases are mixed, 
it turns out that this is equivalent to a positive-only mix- 
ture but with a different expansion parameter fi'^ rather 
than fia, i-e. the function is equivalent to expanding in 
terms of positive powers of (fi — jJ-'a) in which case we 
again get the tracking behavior discussed above. 

A useful way to visualise the results of this paper is 
through the two phase plots in Figsl?! and [SI For each 
plot we pick a range initial conditions and we look at the 
function ^^^(/i) for the choice of fig such that fi^ = 0.4 
in the tracker limit. In Fig. 5 (Bekenstein's model) we see 
that all the curves join the tracker solution. In Figl5]we 
plot a model with n = 2 and n = —1. With some choices 
of initial conditions we have a temporary tracker due to 
the n = 2 factor and eventually rolling to the singularity 
mentioned in the section on the n = — 1 case. While with 
other choices of initial conditions, the tracker is avoided 
all together. 

The next step in the analysis of the cosmology of TeVeS 



is to proceed with the programme initiated in 's'l . In that 
work, it was shown that, for Bekenstein's choice of func- 
tion that it was possible to flt observations of large scale 




a 3 4 5 



FIG. 8: Phase plot for the n — 2 and n = — 1 model. The 
upper curves don't join the tracker solution {Q^=0.4) and 
tracking becomes unstable eventually leading to the singular- 
ity mentioned in the text (in this case with a (f) dominated 
Universe). We used C2C_i < 0. 



structure and fluctuations in the cosmic microwave back- 
ground (CMB), the later if massive neutrinos of mass 
2eV are present. This neutrino mass is on the bor- 
der of current neutrino mass experiments but not 
yet ruled out. It has been claimed that these theories 
predict low odd peaks in the angular power spectrum of 
the CMB due to the absence of dark matter |ll, [l9|. 
This is incorrect: as explained in Q (and subsequently 
reconfirmed in @) the extra fields in TeVeS sustain the 
growth of perturbations through recombination, driving 
the acoustic oscillations in a manner entirely analogous 
to that of dark matter. Whether the power spectra based 
on a generic function along the lines considered in this 
work can fit current observations is left for a future work. 
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